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ABSTRACT 



Building a self-gravitating hydrodynamic code as a combination of a hydrodynamic 

^ . solver and a gravity solver is discussed. We show that straightforward combining those two 

0\ . solvers generally leads to a code that does not conserve energy locally, and instead a special 

^ ! Consistency Condition ought to be satisfied. A particular example of combining Softened 

, Lagrangian Hydrodynamics (SLH) with a P^M gravity solver is used to demonstrate the 

effect of the Consistency Condition for a self-gravitating hydrodynamic code. The need to 

^ ! supplement the SLH method with the P^M gravity solver arose because the Moving Mesh 

! Gravity solver, used in conjunction with the SLH method previously, was found to produce 

^ ! inaccurated results. We also show that most existing cosmological hydrodynamic codes 

[ implicitly satisfy the Consistency Condition. 
^ . 
O ■ 
(N ; 

I Subject headings: cosmology: large-scale structure - cosmology: theory - dark matter - 

^ ' hydrodynamics - numerical methods 

^ : 
o ■ 

^ , 1. Introduction 

■ 

' Numerical cosmological simulations are now widely used to assess properties of various cosmological 
' models and to study theoretical aspects of large-scale structure evolution. The variety of numerical 
' methods employed ranges from high resolution coUisionless A^-body simulations (e.g. Davis et al. 1985; 
■ ^- ' Park 1990; Bertschinger & Gelb 1991; Gelb & Bertschinger 1994a, 1994b; Xu 1995) to sophisticated gas 
dynamics methods (e.g. Cen & Ostriker 1992a; Katz, Hernquist, & Weinberg 1992; Evrard, Summers, 
& Davis 1994; Bryan et al. 1994; Summers, Davis, & Evrard 1995) to numerical algorithms that at 
some level include galaxies together with coUisionless dark matter and intergalactic gas as a distinct 
third component of a simulation (Cen & Ostriker 1992b, 1993a, 1993b; Frenk et al. 1995; Gnedin 
1996a, 1996b). 

However, not all of these simulation methods are reliably tested and well understood. Historically, 
coUisionless A^-body methods have been subjected to the most detailed study and testing. The theory 
of A^-body methods is well developed and their errors and limitations are well understood (Hockney & 
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Eastwood 1981; Efstathiou ct al. 1985). Cosmological hydrodynamics methods range from Eulerian 
hydrodynamic techniques borrowed from engineering appUcations (Cen et al. 1990; Ryu et al. 1993; 
Bryan, Norman, & Ostriker 1995) to quasi-Lagrangian methods that include the "Smooth Particle 
Hydrodynamics" (SPH) method (Evrard 1988; Hernquist & Katz 1989) and "Moving Mesh approach" 
(MMA) (Gnedin 1995; Pen 1995b). Unfortunately, currently there is no full understanding of errors 
and specifics of various hydrodynamic methods. The first comparisons between various cosmological 
methods (Kang et al. 1994; see also Gnedin 1995) demonstrated that differences between various 
methods can be so substantial that only a few statistical quantitative results are consistent between 
different methods. Further work on comparison of various numerical techniques is therefore required 
to support quantitative results of cosmological hydrodynamical simulations. 

Since the Moving Mesh approach is the most recent and therefore least understood, it requires 
special attention in comparing it with other existing numerical techniques. In this paper we 
concentrate on comparing an implementation of the Moving Mesh approach called Softened Lagrangian 
Hydrodynamics (SLH; Gnedin 1995) with existing A^-body techniques. One possible advantage of the 
Moving Mesh approach is that it offers a way to calculate gravitational forces with high resolution 
without employing a "particle" approach like P^M or TREE. Since the moving mesh represents 
a single-valued coordinate transformation from quasi-Lagrangian space into real (physical) space, 
the Poisson equation in real space can be reduced to an elliptic partial differential equation in 
quasi-Lagrangian space, which can then be solved by the standard numerical techniques. In the SLH 
method thus calculated the gravity force has also the virtue of having exactly the same resolution as 
the hydrodynamic solver has (we elaborate below what this actually means). There are also efforts 
under way to develop a purely A'"-body version of the Moving Mesh approach (Pen 1995a). 

The Moving Mesh approach looks very promising since it has higher spatial resolution in dense 
regions than Eulerian codes for the same amount of computational resources and is substantially 
faster than SPH methods for the same spatial resolution. We, therefore, undertook to investigate its 
accuracy by comparing SLH with the P^M gravity solver. Our detailed tests however have uncovered 
two serious problems with the Moving Mesh gravity solver which have not been discussed in the 
literature; the problems we found may appear in a wide range of algorithms, so we report them in 
this paper. The first problem is the discreteness noise (shot noise) in high density regions. It is 
well known that the shot noise does not represent a serious problem in traditional gravity solvers 
(Hernquist, Hut, & Makino 1993), where force errors effectively average out in high density regions; 
however, in the Moving mesh approach the relationship between the density and the gravitational force 
is nonlinear, with the nonlinearity coming from the nonlinear dependence of the volume of a mesh 
element on the force acting on it. In this case the shot noise does not necessarily average out as in the 
linear case, leading to spurious numerical heating of high density regions. Secondly, a solution to the 
finite- difference realization of the Poisson equation on a strongly deformed mesh does not necessarily 
converge to the true solution, giving the incorrect force law. We propose a solution to the first problem 
in §2, but can not resolve the second one. We have therefore been forced to abandon the Moving Mesh 
gravity solver and have turned instead to a P^M gravity solver that is known to work correctly and 
accurately. However, while developing the new SLH-P"^M code, we encountered another problem which 
is generic for (and potentially present in) any self-gravitating hydrodynamic code, not only a Moving 
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Mesh code: since the gravity force solver and the hydrodynamic solver may treat spatial gradients 
differently, they are not a priori consistent with each other. This inconsistency may lead to a serious 
local nonconservation of energy even if, globally, energy is well conserved. We discuss this effect in 
detail in §3, where we also derive the "Consistency Condition" that every cosmological hydrodynamic 
code must obey in order to be energy conserving and apply it in §4 to develop a new SLH-P^M code. 
Finally, we repeat some of the Kang et al. (1994) tests for the new code in §5. We summarize our 
conclusions in 56. 



2. Failure of Moving Mesh Gravity 

The Softened Lagrangian Hydrodynamics method was extensively tested in Gnedin (1995); 
however, the main emphasis there was on the hydrodynamic properties of the code. Later, we 
proceeded to test the gravitational solver of the SLH code and discovered interesting problems with 
the Moving Mesh gravity solver and potentially with self-gravitating hydrodynamic solvers in general. 
The purpose of this paper is to report the results of those tests and our solutions to the problems. 

We have chosen a CDM+A model as a framework for our tests. The cosmological parameters we 
use are as follows: Qq = 0.4, h = 0.65, and ag = 0.85. We use the BBKS transfer function (Bardeen 
et al. 1986) to set up initial conditions. The model is COi^ii'-normalized and is in a reasonable 
agreement with all existing large- and intermediate-scale observations (Kofman, Gnedin, & Bahcall 
1993; Peacock & Dodds 1994; Ostriker & Steinhardt 1995). We have run several simulations with 32^ 
and 64^ grids all of them having the initial (Eulerian) comoving cell size (or, equivalently, the mean 
interparticle separation) equal 0.5h~^ Mpc, and, thus, a 32^ simulation has a 16h~^ Mpc box size and 
a 64^ simulation has a 32h~^ Mpc box size. We set both the P^M and SLH softening parameter to 
1/10 (the softening length is 50/i"^kpc) and keep it constant throughout our tests. Our simulations 
have box sizes too small to be of interest for large-scale structure formation at the current epoch, but 
they have sufficient resolution to push simulations into the highly nonlinear regime to test the SLH 
method in challenging conditions. In addition, we assume that the gas is "adiabatic", i.e. no radiative 
processes like cooling or ionization are taken into account. Our current goal is to investigate in detail 
gravitational properties of the SLH code and, therefore, we do not need to complicate our tests with 
radiative physics. 

We mostly concentrate on properties of the largest clumps found in our simulations, unlike Kang 
et al. (1994), who mainly concentrated on some average statistical properties of gas simulations with 
different sizes 0. First, we have run SLH and P'^M 32'^ simulations with identical initial conditions and 
with fib = for the SLH code (the SLH code is thus used as a collisionless iV-body solver, however, 
the hydrodynamic part is used to follow the mesh deformation as the mesh is tied to baryonic gas in 
the SLH method). Fig.|l| shows dark matter particles in a subvolume of these simulations containing 
the largest clump at current epoch. The lower left panel shows particle positions for a SLH simulation, 
and the lower right position shows P^M particles. One can easily see that the SLH method fails to 
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Fig. 1. — Dark matter particles in a subvolume of four 32^ simulations: the original SLH Qf, = (dark matter 
only) simulation (lower left panel), the pure P^M simulation (lower right panel), the original Qx = (baryons 
only) simulation (upper left panel), and the fi;, = SLH simulations with dark matter density smoothed in 
quasi-Lagrangian space (upper right panel). 

produce a clump nearly as dense as the P^M clump. In our comparison we explicitly assume that the 
P^M gives results that are correct for dark matter (Qi, = 0) up to the softening length and we use it as 
a template against which to compare the SLH results. 

For comparison, we have also run the same SLH simulation but with Qb = ^o, which implies that 
the dark matter does not gravitate. The corresponding subvolume is plotted in Fig.Q, upper left panel, 
and, as can be easily seen, there is much better agreement between the SLH and P^M simulations in 
that case 0. We must therefore ask ourselves what is the difference between the SLH = and SLH 
fix = cases that causes such a significant deviation. 



*Thc SLH gas feels isotropic pressure while the dark matter has an anisotropic pressure tensor, so some differences are 
expected. 
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There are two principal differences between dark; matter only and gas only simulations with SLH. 
First, the dark matter density assignment uses the Cloud-In-Cell method to assign the dark matter 
density on a quasi-Lagrangian mesh (as explained below), and this assignment suffers from discreteness 
noise (or, in other words, shot noise). Second, the baryonic pressure tensor is isotropic, whereas the 
dark matter pressure tensor is not. However, the second reason can not explain why a cluster failed to 
form with = SLH, since in the P^M simulation the dark matter has anisotropic pressure tensor as 
well, but a dense cluster forms. We conclude therefore that it is the SLH density assignment scheme 
that is responsible for the cluster diffusion problem. 

In order to elaborate on this effect, we recall the basic ingredients of the Moving Mesh approach 
in general and the SLH method in particular (the reader should refer to Gnedin (1995) for complete 
explanation of notation and terminology). Let us consider a general coordinate transformation, which 
connects the real space coordinates, x', with some other coordinate system g*^, which we will call 
"quasi-Lagrangian" hereafter, 

x' = x\t,q'). (1) 

There are no restrictions on g*^ up to now, but we will assume in the following that the new coordinate 
system q'^ has the property that it redistributes numerical resolution into high density regions. The 
important quantity that controls the degree of deformation between coordinate systems and q'' is 
the deformation tensor, 

dq^ 

Then the density of gas or dark matter can be represented as 



4 ^ ^- (2) 



p = po/A, (3) 

where A is the determinant of the deformation tensor A^, 

A = detA^„ (4) 

and Pq = Po{q^) is the quasi-Lagrangian density, or, in other words, the mass of a fluid element having 
a unit volume in quasi-Lagrangian space q'^. Note that quasi-Lagrangian coordinates may deviate 
significantly from exactly Lagrangian ones; i.e., for a given mass element q'^ need not be constant. 

In a Moving Mesh approach (and, therefore, in the SLH code as well), coordinates x^ are real space 
images of vertices of a uniform mesh in quasi-Lagrangian space, qjjf. — i, qfj^. — j, qfji^ — k. Since A, as 
a function of the deformation tensor, is the property of the mesh, it is a smooth function of position (as 
smooth as the actual mesh equation allows it to be). For baryonic gas, po,B is also a smooth function of 
position since it represents the mass of the gas in a cell; for a fully Lagrangian flow po,B = 1- However, 
this is not necessarily true for the dark matter density, which may, therefore, have large gradients in 
g-space. In the SLH approach the dark matter density that sources the right hand side of the Poisson 
equation is determined by assigning dark matter mass to cells in the quasi-Lagrangian space with the 
cloud-in-cell technique: 

Po,x{q') = E mi,j,kW{q' - q[j^^). (5) 
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In voids where, by construction, there are always almost the same number of dark matter particles per 
cell (since in voids the baryonic flow is Lagrangian and the pressure is negligible, dark matter velocity 
is equal to the gas velocity at the same position and, therefore, both the dark matter and baryons move 
simultaneously in voids), and, therefore, po,x is a smooth function of position and is almost constant. 
(It would be exactly constant before dark matter shell crossing and before gas shocking if the pressure 
in voids was zero and quasi-Lagrangian space was exactly Lagrangian.) This is not so in high density 
regions. Dark matter particles wander around in a cluster while the mesh does not move substantially. 

In Eulerian PM codes the density fluctuations are a less severe (while still existing) problem 
since there are many dark matter particles per one PM cell in a high density region and the dark 
matter density is determined accurately with the Cloud-In-Cell method. However, in the Moving Mesh 
approach, if the number of dark matter particles is equal to the number of baryonic cells, there is 
on average one (or slightly more than one if the dark matter is more concentrated than the baryons) 
dark matter particle per cell even in the highest density region since the high density is achieved not 
by collecting many particles in one immobile cell but by shrinking a cell to small volume keeping its 
mass approximately constant. If these particles are distributed randomly with respect to the mesh, i.e. 
if all ql j i. are random positions, equation would predict that po,x fluctuates on a scale of one cell 
with the amplitude comparable with its mean value. That means that the density px in clusters would 
fluctuate with the amplitude comparable with its mean value even in the cluster center. Since particles 
move, the dark matter density will fluctuate both spatially and temporally. The temporal fluctuations, 
however, seem to be less serious since the gravitational time scale in clusters is usually larger than the 
hydrodynamic Courant time scale. 

Such enormous fluctuations will produce fluctuations in the gravitational potential; it has been 
shown (Hernquist, Hut, & Makino 1993) that those fluctuations effectively average out for traditional 
gravity solvers. However, in the Moving Mesh approach the situation is different. Since the mesh 
is usually driven by gravity force, and the gravity force itself is calculated on the mesh, there is a 
nonlinear dynamical relation between the gravitational force and, say, the density defined on the 
mesh (with nonlinearity coming from the nonlinear dependence of the deformation tensor on the 
gravitational force). Therefore, the traditional analysis (Hernquist, Hut, Sz Makino 1993), where the 
relation between the force and the density is linear, does not apply and force errors do not necessarily 
average out as in the linear case. The nonlinearity between the density and the force also serves as 
a dissipative numerical mechanism, which leads to eventual blowing up of a cluster. We would like 
to stress here that this effect is a direct consequence of the shot noise effect in clusters, when in the 
Moving Mesh approach there are only a few dark matter particles per cell even in the highest density 
regions. 

It is easy to demonstrate that this shot noise effect is indeed present in Moving Mesh gravity 
solvers. First, we noticed above, that this effect is absent when the baryons instead of the dark matter 
is gravitationally dominant, and po.B is always a smooth function of position. Second, we can define 
Po,x as a smooth function of position if we smooth it in quasi-Lagrangian space. This procedure will 
not significantly degrade the resolution since it is A, not po,x, that carries the high resolution. The 
smoothing will not affect voids since po,x is almost constant in voids; it will affect places where the 
dark matter density is physically different from the baryon density, i.e. near shocks. We can hope 
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Fig. 2. — Density profiles for four clusters shown in Fig.l: P^M only {solid line), fif, = SLH {dotted line), 
Ox = SLH {short- dashed line), and = SLH with smoothing {long- dashed line). The radius has units of 
the mean interparticle spacing; the force softening is R = 0.1. 

however that these regions occupy a small fraction of the total volume, and, what is more important, 
they do not dominate the gravitational potential. 

To test the shot noise hypothesis we have run another SLH simulation where we have smoothed 
the quasi-Lagrangian dark matter density po,x entering the right hand side of the Poisson equation 
with a Gaussian filter with width equal to three cell sizes in quasi-Lagrangian space (we would like 
to stress that this is not equivalent to smoothing the dark matter density, since the factor A remains 
unsmoothed). The resulting particle positions are plotted in Fig.||, upper right panel. One can easily 
see that a cluster is now formed. Fig.^ shows density profiles for the clusters in these four cases after 
the cluster was identified with the DENMAX algorithm (Bertschinger & Gelb 1991). The SLH gravity 
with smoothing is perhaps a factor of 1.5 softer than the Plummer law for these simulations f\. 



^We have also performed the same SLH simulation without smoothing but with 8 and 64 times more dark matter 
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Fig. 3. — Dark matter particles in a subvolume of two 64^ simulations: the 17;, = SLH simulation with 
smoothing (left panel) and the pure A^-body simulation (right panel). 

In order to test the Moving Mesh gravity further, we have performed two 64^ simulations with 
the parameters described above: one is a P^M simulation and the other one is a SLH simulation 
with f2b = and dark matter smoothing as explained in the previous paragraph. Particles around 
the heaviest cluster in the two simulations are plotted in Fig.^ in the right and left hand panels 
respectively. One can see that at higher resolution (64^ vs 32^ in Fig.Q) the SLH cluster again diffuses 
outwards. Fig.^ shows the density profile for the cluster. Obviously, the SLH gravity effectively softens 
at a scale a factor of 6 or 7 larger than the formal softening length. Since the effect of shot noise in 
high density regions is eliminated in the SLH simulation, one has to admit the existence of yet another 
problem with the SLH gravity. 

In order to investigate the accuracy of the SLH Poisson solver we calculated point mass 
gravitational forces for ten randomly chosen positions on the SLH mesh taken from the 64^ SLH 
simulation at z = 0. The corresponding force laws are plotted in Fig.|^, left panel. Since points are 
chosen randomly on a mesh, their effective smoothing lengths are different (for example in a void 
gravitational smoothing is larger than the original cell size since a cell in a void has expanded during 
the evolution); in general, however, there is an acceptable agreement between the SLH point mass force 
and l/i?^ law (dotted line in Fig.^. Now, instead of ten randomly chosen points, we picked up two 
points that lay close to the center of the biggest cluster in the SLH simulation. Their point mass force 
laws are plotted in Fig.^ with solid lines. These laws deviate significantly from the law at 5 to 6 
softening lengths even if these cells have sizes of the order of a softening length. 

What causes these errors in gravity in a cluster? In order to test if it is due to the loss of accuracy 



particles; there is a definite improvement as the number of dark matter particles increases, but the rate of improvement 
is very slow, and even with 64 times more dark matter particles than baryonic cells the cluster under consideration is still 
significantly diffused out. 
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Fig. 4. — Density profiles for two simulations shown in Fig. 3: P^M only {solid line) and 
smoothing {dotted line). 



SLH with 



in the finite differencing of the gravitational potential, instead of solving the Poisson equation 

^ 2 ' 

we solved three Poisson equations for each of force component separately, 

3Qoa 85 



A/. 



2 dx^ 



(6) 



(7) 



which eliminates the need to differentiate the potential but requires at least three times the CPU time 
used by the potential solver. The corresponding point mass force law is plotted with the thin dashed 
line in the right panel of Fig.^. While there is some improvement in solving for the force compared to 
solving for the potential, it is obviously insufficient to improve the Moving Mesh gravity solver. We 
specifically stress here that in this test, the specifics of the SLH algorithm come only with the strongly 
deformed mesh; the Poisson solver is a generic feature of any Moving Mesh approach and it is this 
that fails. We, thus, have encountered another problem with the Moving Mesh gravity: when the mesh 
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Fig. 5. — Left panel, Green functions for ten randomly chosen points in the 64^ SLH simulation (solid lines) 
and the inverse square law (thin dotted line). Right panel, Green functions for two points chosen near the center 
of the cluster from the left panel of Fig.3 (solid lines) and the Green function for one of those points calculated 
by solving a Poisson equation for three components of the gravitational force on the deformed mesh (dashed 
line); shown with the thin dotted line is the inverse square law. 

is strongly deformed, the Poisson solver occasionally gives incorrect results. It seems quite difficult to 
find out what specifically is wrong with the Poisson solver or what properties of the mesh cause such a 
big error. It is also not clear if the less deformed mesh has the same problem but less pronounced so 
that it avoids clear recognition. 

Even if the problem we have highlighted appears only in small volume of the whole simulation, 
there is no guarantee that the rest of a simulation is correct or has a specified resolution. Until the 
specific problem with the Moving Mesh gravity is identified and resolved, we feel that it is unsafe to use 
the Moving Mesh gravity solver for cosmological simulations. The full analysis of Poisson solver errors 
on strongly deformed meshes is beyond the frame of this paper and, therefore, we have decided to 
abandon the Moving Mesh gravity solver and proceed into combining the SLH hydrodynamic method 
with the P^M solver which is known to give correct results (Bertschinger 1991). We shall see that 
merging these codes introduces another problem to solve. 
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3. Consistency Condition 

We now imagine building a cosmological hydrodynamic code by combining separate hydrodynamic 
and gravity solvers together. Both methods can calculate hydrodynamic and gravitational forces on a 
resolution element scale quite differently, and there is a priori no guarantee that these two methods 
would be consistent with each other. 

Let us consider an example. Let us imagine that we constructed an isothermal gas sphere; 
we utilize so many resolution elements (cells or particles), that we do not have to worry about 
relaxation processes. We now insert this object into a hydro+gravity code which has gravitational 
force softening much shorter than the pressure force softeningQ The object, which would be in real 
physical equilibrium (i.e. in equilibrium with exact force laws for both gravity and pressure force), 
would collapse further since the softened gravity force would be stronger than the softened pressure 
force due to larger pressure force softening. The object would thus form a small and dense clump that 
is a purely numerical artifact. This process is an example of a more general inconsistency between 
hydro and gravity solvers that we discuss in detail below. 

Let us consider the hydrodynamic equations with self-gravity in a general coordinate system q'^ as 
defined above (at this point we do not require that the coordinate system be close to Lagrangian; in 
particular, it may be fully Eulerian when = g*). Let p be the density, be the velocity, P be the 
pressure, E = pvf /2 + P/(7 — 1) be the total energy density, and fi be the gravitational force per unit 
mass of the fluid with the polytropic index 7. The Euler equations in a coordinate system g'^ read: 

(8a) 

= Pofu (8b) 

Pov'fi, (8c) 
where the deformation tensor has been defined above, Bl is the matrix inverse to A]., 

BjAj. = AjBl = 

A = dei{A\) is the determinant of the deformation tensor, which is equal to the volume in real 
space of a fluid element with the unit volume in quasi-Lagrangian space g'^, d?x = Ad^q, po is the 
quasi-Lagrangian mass density, 

Po = pA, (9) 
Eq is the quasi-Lagrangian total energy density, 

Eo = EA, (10) 

^This should be considered only as a thought experiment; in reality hydrodynamic codes have pressure softening 
comparable to the resolution element size and, thus, our assumption of having many resolution elements per pressure or 
gravity softening length is unrealistic. To implement it in practice, one would require to smooth both pressure and gravity 
force over a scale encompassing many resolution elements. 



^ + £i(pow^') = 0, 



opo ^ d 
dt dq 

dpoVi d 



and 



dt dq 
dEo d 



+ (^p,v,w' + AB^P 



dt dq 
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and 

= B^iv' - x'), (11) 

with the dot standing for the partial time derivative, x^t, g*^) = dx'^{t, q^)/dt. 

If we now consider a hydrodynamic code that solves these equations, we must interpret the 
spatial derivative d/dq^ as a numerical derivative (for example, a finite difference derivative for a 
mesh-based code), and we denote it as d^/dq^^ and we must include numerical dissipation terms 
responsible for shock handling and stability. Let us introduce a numerical density flux correction JF^, a 
numerical momentum flux correction Qf ^ a numerical total energy flux correction and a numerical 
gravitational energy correction A through the following formulas: 

^ + |i(Po-' + ^')=0, (12a) 



9i dq 
and 

^ + 1^ {e,w^ + P + = po^^Vi + A- (12c) 

Among those four, only J^^ and A will be of interest to us. 

We now want to establish the total energy conservation. The total energy can be represented as 
the sum of the total gas energy, 

S = j Eod^q, 
Jn 

and the gravitational energy, 

2 Jn 

where we understand integrals in a numerical sense and denote this by index N, and is a gravitational 
potential, and is determined by the following equation: 

0(g^) = / G{\x'(q') - x'iqDDpoiqDd^i, (13) 

and G{r) is a Green function for a gravity solver. It is important to note here that this function is 
arbitrary, and, therefore, our consideration is equally applicable to a softened gravity law as well as to 
the 1/r^ law. Moreover, we do not need to specify a softening length for the gravity. It will be shown 
below that the problem of inconsistency does not reduce to the problem of different resolution scales 
for gravity and gas and is more generic. 

We now demand (neglecting cosmological expansion for the moment) that 

±{S + W)^0. (14) 

It can be easily seen that 



- 13 - 



provided the numerical integral of a numerical divergence is zero (we assume that this is the case). For 
the time derivative of the gravitational energy, we obtain: 



dt 



N at 



N OX^ 



G{\x^-x\ql)\)po{q\)d\, 



N 



d^q. 



(16) 



It is important to note that the derivative with respect to x* that enters the last equation is a full 
derivative , not the numerical one. We can introduce the new quantity, /j, defined as: 



d(j) 

dx"^ 



d 



dx^ 



Gi\x'-x\qf)\)poiq'l)d\ 



N 



d\, 



(17) 



which is the true gradient (with the minus sign) of a numerical gravitational potential, i.e. the gradient 
which will be computed in a simulation, in which all spatial gradients are computed mathematically 
exactly (for example, a P^M code possesses this property). Note, that due to the finite number of 
resolution elements in a real simulation, /j is not necessarily equal to /j, but rather is a numerical 
approximation to it. Equation (|T6|) can now be written as follows: 



d 

—W 
dt 



^P^ A. ^3 

^4>d q 
N at 



N 



pox'fid^q. 



(18) 



We now assume that the numerical derivative obeys the product rule, i.e. 

dq^ ~ ^ dq^ ^ dq^ ' 



(19) 



This is true for dimensionally split mesh codes with central differences and can be true in SPH 
codes with appropriate choice of weightening scheme; it is generally not true in the SLH code since 
it is not dimensionally split, and, therefore, the SLH code has additional energy errors compared 
to dimensionally split mesh codes. (These errors are typically small and random and seem not to 
contribute significantly to the total energy error, although, in principle, nothing prevents them from 
causing a significant energy error.) 



Then we can use equation (|12a|) to reduce equation (|T8|) to the following expression: 



d 

di' 



{pow' + T')^d'q- I pox^f.d'q. 



(20) 



Combining both equation (|T5|) and equation (pO]) with equation (p!^, we finally obtain: 



N \ oq'^ 



) d^q + 


/ Po 




Jn 



v'h-x'U + {v'-x')Bt 



dq^ 



d-'q = 0. 



(21) 



We can immediately derive the form for A from the equation (^) taking into account that A is 
the numerical correction to the gravitational energy change and it must vanish when JF^ vanishes: 



A 



dq'^ 



(22) 
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Now we are left with the following equation, which must be satisfied if the energy is to be conserved: 



Po 

N 



d^q = 0. (23) 



Let us examine this carefully. There are three gravitational forces that enter this equation: the original 
force that a full gravity solver provides, /j, the force that the hydrodynamic code feels, fi, and the 
numerical gradient of the gravitational potential, B^dN4>/dq'^ . The last expression is the gradient of a 
scalar function as the hydrodynamic code understands it; for a mesh-based code it is a finite difference, 
for a SPH-like code it is the function values weighted by the gradients of the kernel, etc. 

If we require that the self-gravitating hydrodynamic code should conserve energy everywhere and 
for all possible cases, we must satisfy the Gravitational Consistency Condition at every resolution 
element: 

v^f. = - K - i')B^^. (24) 



Strictly speaking, we may add a divergence term, which always can be hidden into Ti.^. This condition 
should be considered as a restriction that the gravitational force fi should satisfy in order for the 
self-gravitating hydrodynamic code to be energy conserving. 

It is instructive to consider some special cases: 

1. A fully Lagrangian code (in particular, SPH). For a fully Lagrangian code = and equation 
( plD is satisfied if fi = fi. Thus, any gravitational solver is consistent with the fully Lagrangian 
code; in other words, a fully Lagrangian code is a priori consistent with a gravity solver and no 
special care should be taken to satisfy the Consistency Condition. 

2. A fixed mesh finite-difference code (in particular, a fully Eulerian finite-difference code). For this 
case, we can assume = and the Consistency Condition is satisfied if 

Ji = -^i 



dq^ 

We conclude, therefore, that it is essential that the gravitational force in a fixed mesh code be 
a finite difference gradient of the gravitational potential. In particular, in a PM + Eulerian 
finite-difference code combination it is the gravitational potential that should be calculated in 
the PM part, and not the force. 

3. A Riemann solver. In this case the numerical gradient dj^/dq^ is actually a Riemann solver; it is 
not straightforward to represent the gravitational force as the gradient of a potential in this case; 
rather it is easier to incorporate the gravitational force into the Riemann solver itself as: 



fi = fi = -B- 



dq^ 



where notation di^/dq^ means the Riemann solver with the gravitational force incorporated in it; 
in particular, this is the way the KRONOS code works (Bryan et al. (1995). 
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4. A Moving Mesh gravity code (i.e. a code that uses the Moving Mesh gravity as a gravity solver). 
Strictly speaking, for this case equation (|T3p does not hold, but we can still consider it here 
recalling that 

f - pkdN(p 

for a Moving Mesh Gravity code. In that case fi is also equal to both and —BfdN(p/dq^ and 
the Consistency Condition holds. 

Thus, the Consistency Condition is usually satisfied for most of currently used numerical 
techniques. 

Let us consider here a simple example. Let us imagine that we have a one- dimensional Eulerian 
finite-difference hydrodynamic code (so that the product rule [|19[ is always satisfied) that keeps all 



hydrodynamic quantities at the centers of equal cells numbered by index i = 1,2, N (where N is the 
number of cells and the periodic boundary condition is assumed), so that the gas density at the cell i 
is Pi, the energy density is Ei etc. This code is combined with the PM gravity solver that uses the gas 
density pi to solve the Poisson equation (using, say, the FFT technique) 



+ - 20i = 47rGpi (25) 



(we assume here Ax = 1 for simplicity) to obtain the gravitational potential (pi at the cell centers, 
derives the gravitational force fi acting at the cell center by finite differencing the potential, 

/. = {<P^-l - <f>^+l)/'2, (26) 

and uses it to update the gas velocity t>j. The code constructed this way conserves total energy 
exactly (provided the hydrodynamic part conserves the kinetic plus thermal energy exactly, as is the 
case for conservative schemes) as explained above, i.e. it satisfies the Consistency Condition. Now 
let us imagine that one decides to replace the PM gravity solver with the exact gravity solver, which 
calculates the exact gravitational potential and force / given the density distribution: 



p{x) 



Pi, < a; < 1, 

P2, l<x <2, 

Pi, i — 1 < X < i, 

Pn, N - 1 < X < N 



(One can achieve this, for example, by using the PM gravity solver with a much finer mesh.) If one 
now uses the exact gravitational force fi to update the gas velocity Vi, the Consistency Condition will 
not be satisfied, and local energy errors will be introduced. However, if one uses the exact gravitational 
potential and calculates the gravitational force / using (^) with the new potential, the energy 
conservation will be restored even if the exact gravitational potential is different from the potential 
obtained by solving the Poisson equation on the hydrodynamic mesh. Thus, the Consistency Condition 
requires a consistent treatment of the gravitational potential and forces in the dark matter and gas. 
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We would like to stress here again that the Consistency Condition is not equivalent to having 
gravity and pressure force softening commensurate; the Green function in (|T3p allows for any value 
for the gravitational softening parameter; more that that, since the shape of the Green function is 
not specified, it is not required that the force fi be actually gravitational force; any potential force 
solver (for example, molecular forces) should satisfy the Consistency Condition similar to equation 
(plf). This comment would imply that the requirement of having gravitational and pressure softenings 
commensurate is an additional requirement a self-gravitational hydrodynamic code must satisfy. 



4. A SLH-P^M Code 

We now turn to constructing a SLH-P^M code, where the gravity force is calculated by a P^M 
solver and the mesh velocity is not necessarily equal to the fluid velocity f*. In this case it is 
impossible to solve equation (^) exactly, and only an approximate solution can be proposed. 

We start by noting that the SLH mesh equation, i.e. the equation that connects and w * in the 
SLH approach. 



(9 d 



4) 



dq^dq^ dq 

is a Galilean invariant modification of the following equation 



(27) 



= {6} - (28) 

(here a* is the softening tensor that vanishes in the Lagrangian limit and approaches Sj in the Eulerian 
limit; for full definitions see Gnedin (1995) or Appendix A). If we consider from equation (^8]) as an 
approximation to derived from equation (|27|) and substitute ( pSf ) into (0), we obtain the following 
(Galilean invariant) equation: 



/. = /.- oi { h + 4§ 1 ■ (29) 



This equation will serve as our definition for the gravity force on the gas in the SLH-P^M code. Both 
fi and are obtained by means of a P^M code and ji is the gravitational force applied to the baryonic 
gas. The forces on the dark matter particles are given instead by ji. Tests also show that the following 
simplified version of equation (|29|) is a good choice as well: 



where a is the maximum eigenvalue of the softening tensor cr/. One can note that, again, in the fully 
Lagrangian {a = 0) case fi = fi, while in the fully Eulerian (a = 1) regime fi = —BfdN(f)/dq^. 



Equations (^) and (|30|) have a remarkable feature that they guarantee the commensurate 
softenings for both gravity and pressure for the SLH code, since the SLH code switches from Lagrangian 
to Eulerian description as density increases and, therefore, the gravitational force at the resolution 
limit is the gradient of scalar function in the way the hydrodynamic solver understands a gradient; 
gravitational and pressure resolutions for the gas are again exactly equal as with the original SLH code. 
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Direct force 




Fig. 6. — Baryonic particles (cell centers) for two 32^ SLH-P^M simulations: without Consistency Condition 
{left panel) and with Consistency Condition [right panel). 

However, since the dark matter particles are acted upon by the force fi, the gravitational resolution 
for the dark matter component need not be (but, of course, may be set) equal to the gravitational 
resolution for the gas component. 

Now we turn to the details of implementing SLH-P^M. There exist two distinct ways to combine 
SLH hydro and P^M gravity codes (and many intermediate variants). 

The first one, which we call a "maximally coupled" code, imphes that dark matter "lives" in 
quasi-Lagrangian space like in the original SLH code; then the P^M part is used to calculate forces 
acting on vertices of the mesh, where the mass in a vertex is found by linear interpolation from the 
cell centers and the mass of a cell is determined by the total mass of dark matter, gas, and, possibly, 
galaxies in the cell. This approach requires calculating the P^M force only on Nb mesh vertices, 
however, since the mesh is tied to baryons, the maximally coupled version fails to follow accurately 
places where the dark matter density differs significantly from the gas density, like cores of big clusters. 

The second variant, which we call a "minimally coupled" code, allows dark matter to evolve in real 
space similar to existing SPH approaches. The mesh is then considered to represent the baryonic gas 
only, and only the gas mass is included in computing the cell mass. This variant of the SLH-P'^M code 
does not suffer from the limitations of the minimally coupled mode, but requires substantially more 
CPU time since computing of the gravitational force acting on both Nb grid cells and A^s: dark matter 
particles is required. For the most common case Nx — Nb this would require twice the PM time 
and four times the PP time of the maximally coupled code. When combined with the hydrodynamic 
part, this averages out to a minimally coupled mode being approximately 2.5 times slower than the 
maximally coupled mode. 

We now proceed to testing the SLH-P^M code. Since, in addition to the mesh softening of SLH, 
the P^M gravity solver introduces a Plummer softening, the new SLH-P^M code has two softening 
lengths. In what follows we set these two softenings equal 1/10 so that our simulation resolution will 
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Fig. 7. — Density profiles for two largest clusters in the 64^ SLH-P^M Ox = simulation for dark matter 
{dashed line) and baryonic gas {dotted line); shown with the solid line is the density profile for the same cluster 
from the pure A^-body (P^M only) simulation with the same initial conditions. 

not exceed that of a pure A^-body template run; in principle, the gravitational force acting on the gas 
can be computed with vanishing Plummer softening since the mesh softening will prevent exceedingly 
small scales from being resolved. Therefore, in the maximally coupled mode, when the dark matter, 
being placed on the baryonic mesh, is subject to the mesh softening, the Plummer softening can be 
assumed arbitrary small thus squeezing the last drop of resolution from the simulation. One must note, 
however, that in that case the resolution becomes even more anisotropic and nonuniform across the 
simulation volume. 

First, we would like to demonstrate the importance of the Consistency Condition. Figure ^ shows 
the gas distribution of the biggest clusters in two 32^ baryons-only {Qx = 0) simulations (cf. Fig.|l|): 
the left panel shows the results for the simulation where no Consistency Condition was applied, and 
the gravitational force acting on gas was directly provided by the P^M part of the code; the right 
panel shows the same cluster in the simulation where the Consistency Condition of the form (|30D was 
incorporated in the code (we remind the reader that the corresponding A^-body result is plotted in the 
lower right panel of Fig.^). It is easily seen that without the consistency condition most of the gas 
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Fig. 8. — Total energy error for the SLH-P^M64 run as a function of scale factor a; the total error fluctuates 
around zero with amplitude 1% (except for a < 0.2 where the energy error is not computed accurately enough). 

sinks to the very center of the cluster thus violating energy conservation. This effect of gas sinking 
toward the cluster center can be solved without Consistency Condition by simply smoothing the 
gravitational force on a scale of a cell size; in that case the gravity force would have approximately the 
same softening as the pressure force. This, nevertheless, would not solve the problem since the cluster 
profile would still be significantly different from the iV-body results due to local nonconservation of 
energy. 

Since the original SLH method failed substantially only at the resolution of a 64^ mesh, we ran a 
64^ Qx = minimally-coupled SLH-P^M simulation (the case Qb = coincides with the pure A^-body 
result). The profiles for the two largest clusters in the simulation shown in Fig.|^ together with the 
pure iV-body result (the solid line). The dashed line shows the dark matter profile and the dotted line 
shows the gas profile. While the dark matter is fairly consistent with the pure A^-body profiles, the gas 
distribution has much larger cores. However, we can not attribute this purely to lack of resolution: 
recent high resolution simulations have also found this effect (Navarro, Frenk & White 1995)[[ We thus 
conclude that the agreement is acceptable and proceed to further test the new SLH-P'^M code. 



5. Comparison with Other Codes 

In this section we repeat all the major tests of Kang et al. (1994) to assess differences between 
the repaired SLH code and other hydrodynamic codes. We performed three simulation with 32^, 64^, 



^In real clusters gas cores are also much larger than the dark matter cores as shown by comparison of X-ray pictures 
and lensing reconstruction (cf. Miralda-Escude 1995); we do not use this fact as an argument since it is not proved that 
our tests possess adequate resolution and physics to simulate real clusters. 
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Fig. 9. — Comparison of integrated quantities: volume averaged temperature (T) (upper left panel), mass 
averaged temperature {T)p (upper right panel), X-ray luminosity (lower left panel, in arbitrary units), and 
rms density fluctuations a (lower right panel), for the SLH-P^M approach (solid circles connected with the 
bold solid line) to the other existing codes: original SLH (stars connected with the bold solid line), CO J (filled 
squares), TVD (small filled circles), PSPH (open circles) and TSPH (open triangles). 

and 128^ mesh sizes with initial conditions of Kang et al. (1994) 0. The total energy error in terms 
of the quantity R defined in Ryu et al. (1993) for the SLH-P^M64 run is plotted in Fig.|. One can 
see that the total energy error fluctuates around 1%, which is better than most existing cosmological 
hydrodynamic codes (see Table 1 of Kang et al. 1994 for comparison). 



^These are: Standard CDM model, fib = 1, 6Ah~^ Mpc box size; initial conditions were defined on a 64'^ mesh; for the 
32'^ simulation every other value along each dimension was taken, and for 128^ simulation linear interpolation was used 
to determine initial values on a 128'^ mesh. Note, that the latter procedure in fact adds some small scale power to the 
initial conditions and it is incorrect to say that the 128^ SLH-P'^M simulation and the 128^ and 256^ Eulerian simulations 
of Kang et al. (1994) have the same initial conditions as the 64'^ simulations. In particular, the lack of convergence in 
the Kang et al. (1994) Eulerian simulations is, due at least partially, to this effect. However, we show results of a 128^^ 
simulation as a comparison with larger Eulerian simulations since the initial condition for our 128'^ simulation are similar 
to 128^ Eulerian simulations of Kang et al. (1994). 
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versions only. 



Fig. 10.— A slice of original (256^ rebinned) data for SLH-p3M64 (a), and SLH-P^M 128(6 ) runs. The top 
row shows density and the bottom row shows temperature; the right column shows a zoom of a 1/16 of the 
whole slice (24 < x < 40, 40 < y < 56) containing one of the largest clusters present in the slice. 

The final data for each simulation were rebinned onto a uniform 16^ mesh, and integral properties 
of the resulting 16^ solution were calculated. Figure ^ compares these against each other and against 
other codes shown in earlier work (Kang et al. 1994; Gnedin 1995). Plotted on four panels are the 
volume averaged temperature (T), the mass averaged temperature {T)p = {pT)/{p), the rms density 
fluctuation a = (p^)/ (p)^ — 1, and the total X-ray luminosity = p^T^^'^ in some arbitrary units. 
Four other cosmological codes are presented together with the SLH code described in this paper; 
two Eulerian hydrodynamic codes: the Cen-Ostriker- Jameson code (Cen et al. 1990; filled squares), 
the TVD code (Ryu et al. 1993; filled circles), and two SPH codes: Tree-SPH of Hernquist & Katz 
(1989; open triangles) and P^M-SPH of Evrard (1988; open circles). The self-evident notation used 
to designate these codes is fully explained in Kang et al. (1994). The bold solid line with open stars 
tracks the results of the original SLH code (which we shall hereafter call SLH-MMG as an abbreviation 
to SLH-Moving Mesh Gravity) and the bold solid line with solid circles shows the results of the new 
SLH-P^M code. One can easily see that the SLH-P'^M code provides a major improvement in resolution 
over the original SLH-MMG code, with resolution similar to or exceeding that of the SPH codes for a 
64^ grid of particles while the temperature is similar to Eulerian codes. The SLH-P^M code, therefore, 
manages to reach the high resolution of SPH codes while introducing less numerical heating. The 
SLH-P^M code is also a factor of 2 to 3 faster than a SPH code for the same resolution. 

Let us now compare original (unrebinned) data for different codes. Fig.[lO|a shows a one-cell-thick 
slice of original data of the SLH-P^M64 run rebinned onto a 256^ mesh (to facilitate the comparison 
with the TVD 256 run); the slice has 256^ points and is 0.25/i~^Mpc thick. Fig.p!0|b shows the same 



slice for SLH-P^M128. Also shown are subslices containing one of the largest clusters found in the 
whole slice. 

These slices are to be compared with corresponding plots from Gnedin (1995), however, in order 
to ease comparison, we plot subslices shown on the right side of FigJTo|a,b together with corresponding 



subslices for TVD 256, PSPH64, and SLH-MMG 64 simulations separately in Fig.|Il] and Fig.|l^; 
both figures show PSPH 64 and TVD 256 results on the bottom; Fig.|ll] shows comparison between 
the original SLH-MMG 64 simulation (upper left panel) and the SLH-P'^M simulation (upper right 
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Fig. 11. — Zooms of the density contours of 1/16 of the whole 256 x 256 sHce (16 < x < 32, 16 < y < 32) 
for four different codes: SLH-MMG,64 {upper left panel), PSPH64 (lower left panel), SLH-P^M64 (upper right 
panel), and TVD 256 (lower right panel): (a) the density distribution; (b) the temperature distribution. 



This figure is too large for ASTRO-PH. The full paper including color 
gif versions of figures 10-12 can be found at 
l:tp:/ / arcturus.mit.edu/ Preprints/ sihp^m.ttp. tar. g^ ; 

the file 
sIhpSm. complement. tar. gz 
in the same directory contains figlO-12 in b/w postscript and color gif 

versions only. 



Fig. 12.— The same as Fig.ll except SLH-MMG,64 and SLH-p3M64 are replaced by SLH-p3M64 and 
SLH-P^M 128 respectively. 

panel) and Fig.|l^ plots SLH-P^M64 and SLH-P^M128 simulations in the upper row to demonstrate 
improvement and convergence. 

One first notes that there are no significant differences visible between the original SLH-MMG 64 
and SLH-P^M 64 results; the SLH-P^M 64 density distribution resembles that of TVD 256 slightly 
better, but shocks seem to be more blurry since the higher resolution of the SLH-P'^M code leads to 
fewer cells being left in voids and in the vicinity of shocks. The comparison between SLH-P^M 64 and 
SLH-P^M 128 shows little difference in the density distribution but substantial improvement in the 
shock thickness. Also the SLH-P^M 128 result has somewhat higher temperatures and shock extent 
because the gas was heated to higher entropy by the extra small-scale power. 

It is quite instructive to compare the SLH-P^M results to SPH and TVD slices. One can notice 
here that, except in the highest density regions, the real resolution of an SPH simulation is low; it is 
hardly possible to identify filaments at all, whereas they are easily visible in the SLH and TVD plots. 
In order to understand this somewhat unexpected conclusion one must recall that the hydrodynamic 
quantities in the SPH method are computed by averaging over, say, 30 nearest neighbors, while in the 
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SLH approach they defined in every cell and do not require any averaging to compute. Since in the 
vicinity of filaments the SLH mesh closely follows the gas fiow, and, therefore, the hydrodynamical 
description is close to Lagrangian, the SPH description of the gas is similar to the SLH description 
being smoothed over 30 neighboring cells, which roughly corresponds to 1/5 — 1/8 of the whole slice 
size depending on the density. More than that, only at overdensities above 30 does the SPH method 
acquire higher resolution than the Eulerian hydrodynamic code with the same number of cells, and the 
TVD 256 run becomes inferior in spatial resolution to the SPH 64 run only at overdensities in excess of 
30 X 64 f« 2000. 

In general, the SLH-P^M 128 run compares better to the TVD 256 run than both SLH-P^M 64 and 
PSPH 64. In some places the correspondence between the former two runs is quite good, but also there 
are places where the two codes differ substantially. We therefore conclude that in detail, differences 
between existing hydrodynamic codes are still substantial, and much work is required in order to 
achieve the same level of computational accuracy and consistency as is observed within coUisionless 
A^-body codes. 

6. Summciry 

We have demonstrated that Moving Mesh gravity solvers suffer from two serious limitations: the 
shot-noise in high density regions and incorrect Green functions. We have proposed a way to solve 
the former problem and failed to completely identify the source of the latter. Since our Moving Mesh 
gravity solver can not be trusted now, we have proceeded further by building a combined SLH-P^M 
code, where the SLH hydrodynamic solver is combined with the P^M method for computing gravity. 

In order for a self-gravitating hydrodynamic code to conserve energy, a special "Consistency 
Condition" ought to be satisfied. We have derived the Consistency Condition for the SLH-P^M code 
and showed what errors would appear should this condition be violated. Ignoring this condition can 
lead to serious errors in dense regions. 

Finally, we have compared the SLH-P^M code with existing cosmological hydrodynamic codes 
including the original SLH code. The SLH-P^M approach offers a substantial increase in resolution, 
making it comparable to the SPH codes in the highest density regions, while maintaining reasonable 
accuracy in resolving shocks and hydrodynamic features. In particular, we have repaired the biggest 
problem of the original SLH method: failure to simulate clusters accurately, which stem out of a faulty 
Moving Mesh gravity solver. 

We, therefore, conclude that the new SLH-P'^M code is a good numerical tool to investigate the 
formation of nonlinear structure in the universe including both gas and dark matter. 

We are thankful to J. P. Ostriker for many fruitful discussions; parallelization of the SLH-P^M code 
for a shared memory multiprocessor was improved with suggestions from F. Summers. Supercomputer 
time was provided by the National Center for Supercomputing Applications. This work was supported 
by NSF grant AST-9318185 awarded to the Grand Challenge Cosmology Consortium. 



-24- 



A. Softening Tensor 

At a general point in the flow the softening tensor is a function of the deformation tensor A], fj^ let 

g^^ ^ Al6'^Ai (Al) 

be the conjugate metric of quasi-Lagrangian space and let \g{jj^^) be its eigenvalues. Since g^^ is a real 
symmetric matrix, one can introduce a unitary matrix C* such that 

-g^^ = Cld\^g[Xgf CI (A2) 

where diag [A^]"^ denotes a diagonal matrix with entries. 

We can now introduce the softening tensor a*-' built according to the following formula: 

r I -I ab 

a'^ = C:di^g[cri^g)] CI, (A3) 
where cr(A) is defined by the following equation: 

where A* is a softening parameter. In two opposite limits, when all three A^ approach zero or infinity, 
the softening tensor tends to zero or 6^^ correspondingly. 
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